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Abstract 

The contribution of the a 3 single electron-loop vacuum-polarization diagrams 
to the Lamb shift of the muonic hydrogen has been evaluated recently by two 
independent methods. One uses the exact parametric representation of the 
vacuum-polarization function while the other relies on the Pade approxima- 
tion method. High precision of these values offers an opportunity to exam- 
ine the reliability of the Monte-Carlo integration as well as that of the Pade 
method. Our examination covers both muonic hydrogen atom and muon g — 2. 
We tested them further for the cases involving two-loop vacuum polarization, 
where an exact analytic result is known. Our analysis justifies the result for 
the Lamb shift of the muonic hydrogen and also resolves the long-standing 
discrepancy between two previous evaluations of the muon g — 2 value. 
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I. INTRODUCTION 



In a recent paper we evaluated the contribution AE^ of the sixth-order electron-loop 
vacuum-polarization diagrams to the muonic hydrogen ([i~p + ) Lamb Shift [|TJ. Together 
with the proposed measurement of the 2Py 2 — 2Sy2 Lamb Shift || it will lead to a very 
precise determination of the proton charge radius. 

The most laborious part of this calculation is that of the single electron loop diagrams 
contributing to the sixth-order vacuum-polarization function (see Fig. 1). Using the para- 
metric representation for this function || we find its contribution to the Lamb shift to be 

AE ip6) = 0.017 410 (9) m r {Zaf{^j\ (1) 

where m r is the reduced mass of the muon-proton system and Z — 1 is the proton charge 
in units of the electron charge |e|. As a cross-check, we have also computed AE^ using 
the Pade- approximation of the single electron loop sixth-order vacuum-polarization function 
n 3 (z), z = q 2 /4m 2 ., where m e is the electron mass Inserting the real part of LT^ derived 
from Eq. (7) of Ref. [§J into Eq. (3) of Ref. we obtained the Lamb Shift contribution 
in the [3/2] and [2/3] Pade approximations. Since they are practically indistinguishable, we 
will not label them separately and simply quote them as 

AE%? = 0.017 414 9 (25) m r (Za) 2 ^ 3 . (2) 

We also derived the imaginary part of the approximate LT^ and, inserting it in Eq. (6) of 
Ref. fill, obtained 



AE\ P ^ = 0.017 414 9 (26) m r (Za) 2 [ — ) . (3) 



a 
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Again the [3/2] and [2/3] Pade approximations are nearly indistinguishable. The results (0) 
and (|3]) are consistent with each other and agree within one standard deviation with the 
direct calculation (H). Note, however, that the uncertainties in (H) and (|3|) are those caused 
by numerical treatment of the Pade approximation and do not represent the accuracy of 
the Pade method itself. To gain insight in how good the Pade approximation is we have 
examined two cases where exact results are known. Based on these results we argue that 
the true value of the muonic hydrogen Lamb shift will be found within 0.000 7 percent of 
the Pade value, which is well within the uncertainties quoted in (0) or (Kl). 



II. DERIVATION OF PADE APPROXIMATION RESULTS 



Let us begin with a brief review of the derivation of the Pade approximation to the sixth- 
order vacuum-polarization function LT^ (2;), z = g 2 /4m 2 , by Baikov and Broadhurst jJJ. The 
analytic properties of II3 (z) they utilized are: (a) The first three coefficients of the Taylor 
expansion around z = 0. (b) The first two log ^-dependent coefficients of the expansion 
in 1/z for large negative z. (c) The threshold Coulomb behavior which is determined by 
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FIG. 1. Sixth-order vacuum-polarization diagrams with a single electron loop. 



nonrelativistic quantum mechanics. Taking account of these informations they constructed 
the function 



if (2) = ETO + 4n 2 (^) + (1 - z)G(z)[ -G{z) + ^ + 
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16 322 



229 173 



32z 96 



(4) 



where n 2 (z) is the fourth-order vacuum-polarization function ||. G(z) is the hypergeometric 

3. 

2 ! 



function 2^i(l, 1; §; z). For the negative real z, G is given by 

1 



G{z) 



ln(v^i + a/1 - z). 



z 2 — z 

Analytic continuation of G(z) from z < to z > 1 through the upper z-plane yields 

-1 



ReG(z) 
ImG(z) 



\J z 1 — z 
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ln(v^+ y/z- 1), 



2^ 



z* — z 



(5) 

(6) 
(7) 



By construction 63^(2:) is analytic in the z-plane cut along the real axis 1 < z < 00. The 
function defined by 



l-u 



M 1] (-oc 



4u 



1 + w 



(8) 



is analytic for \u\ < 1, and may be simulated accurately by a Pade approximant P{oj). 
Analytic informations on nj^ listed above are translated into six data for P(o>): {P(— 1), 
P(0), P'(0), P"(0), P'"(0), P(l)}. Using these data we have constructed 



a + aidJ + a2cu 2 + a 3 u 3 



(9) 



The coefficients a's and 6's for both [2/3] and [3/2] Pade approximations calculated from 
the given data are listed in Table I. Once a Pade approximant is constructed, we can readily 
obtain the corresponding IT^ from (^) and (0). 

The 2P1/2 — 25*1/2 Lamb shift of the muonic hydrogen has been evaluated in two ways. 
One uses the formula which contains IL(z) for negative z [jlj: 
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TABLE I. Coefficients of Pade approximants. We set 6q = 1 for the overall normalization. 



Coefficient 


[2/3] Pade 


[3/ 2] Pade 


On 


5.450 103 092 


5.450 103 092 


CLi 


-0.966 458 776 


-0.891 171 812 




-1.785 150 929 


-1.800 086 980 


a 3 





-0.025 240 917 


bo 


1 


1 


bi 


-0.456 709 419 


-0.442 895 559 


b 2 


-0.121 731 656 


-0.128 331 492 


h 


0.001 706 927 






where 



and 



2 r°° 

AE = -{Zafm r / ckp(a 2 )n(-a 2 /(4/? 2 )), (10) 

7T JO 



2x 2a 2 (l - a 2 ) 



777 

(3 = 1- = 0.737 383 76 (30) . (12) 

m r Za 



Substituting II3 (z) determined from Table I in Eq. ( |i~0| ) and evaluating it numerically [], 
we obtained the result (0). 

The second approach utilizes the imaginary part of the approximate II^ for z > 1 
obtained from the Pade approximant P{ui) by taking its value on the unit circle in the 
upper-half u plane. Using this information the 2P\/2 — 2S\/i Lamb shift of the muonic 
hydrogen can be expressed by |] 

AE = mJZa) / dz . — . A . 13 

V 1 h 7T (1 + 2/V^) 4 V ; 

The result (|3|) follows from fli~3|). 

III. DISCUSSION 

The results (|2|) and (|3]) obtained by the Pade approximation method are in good agree- 
ment with the direct result ([[]). The difference between them is within one standard deviation 



1 Unless specified otherwise integrals are evaluated numerically on DECa using the adaptive- 
iterative Monte-Carlo subroutine VEGAS W\. 
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of the result ([J) and can be ignored for the purpose of comparison with experiment. How- 
ever, (fj), (g) and (|3|) all involve some uncertainties inherent to their derivation. Thus it will 
be worthwhile to examine the nature of these uncertainties. There are at least two possible 
causes which may contribute to these uncertainties: One is that the error estimate generated 
by VEGAS in the evaluation of ([1]) might be a gross underestimate of the true error. The 
other arises from the fact that the Pade approximant does not represent nj^z) accurately 
for all values of z. 

A. Non-statistical error in VEGAS calculation 

The integration routine VEGAS is an adaptive-iterative procedure based on random 
sampling of the integrand ||. In the z-th iteration, the integral is evaluated by sampling 
the integrands at points chosen randomly according to a distribution function chosen in the 
(i — l)-st iteration. This generates the approximate value /j of the integral, an estimate of its 
uncertainty <7j, and the distribution function to be used for the next iteration. After several 
iterations li and Oi are combined under the assumption that all iterations are statistically 
independent. The combined value and error are given by 

^ = (E(^M 2 ))/(E( 1 M 2 )), 

- = (E(iM 2 )r 1/2 * (14) 

i 

In general VEGAS is found to converge rapidly for a sufficiently large number of sam- 
plings. However, a special care is required when it is applied to the integration of Feynman 
amplitudes in which renormalization is carried out on the computer relying on point-by-point 
cancellation of singularities between the unrenormalized integrand and the corresponding 
renormalization term. This does not pose a problem if we operate with infinite precision. 
In reality, however, calculation is carried out using finite precision arithmetic, such as dou- 
ble (real*8) or quadruple (real* 16) precision. Random numbers generated by VEGAS for 
sampling of the integrand will inevitably hit points very close to some singularity. This will 
result in evaluation of the difference of two very large and nearly identical numbers with 
finite amount of significant digits. At such a point most of the significant digits cancel out 
leaving only a few significant digits or no significant digit at all. (This will be referred to 
as digit-deficiency or d-d problem, and the subdomain of integration where this happens, 
which contains some boundary surface of the hypercube in our problem, will be called the 
d-d domain.) This introduces nonstatistical noise in the evaluation of the integral and its 
error estimate, even though the effect, being confined to the d-d domain of very small mea- 
sure, is often not readily distinguishable from the fluctuations inherent to random sampling 
of the integrand. They tend to give the integral a false value and might cause deviation 
of error estimate from the assumed statistical behavior. However, the d-d problem encoun- 
tered in the z-th iteration will not affect the performance of the (i + l)-st iteration unless 
it distorts the distribution function very severely. As is seen from (|l4l) J, with larger are 
given smaller weights, giving only small impact on the composite / and a. Problem arises, 
however, if some a,i is relatively small even if it is suffering from the d-d problem. In such a 
case the final / may be distorted in an unpredictable way. 
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Relative impact of the d-d domain decreases as the sampling statistics in an iteration 
increases. More importantly, it decreases dramatically if quadruple precision is adopted. In 
the past, however, this was not necessarily a practical approach because it typically requires 
20 or more computing time for execution compared with double precision calculation. Only 
recently the availability of faster computers has made this a viable option. 

Because of higher speed we normally start evaluation of Feynman integrals in double 
precision. If this runs into a d-d problem, we split the integration domain into a small region 
around the d-d domain and the remainder. The difficult region is then evaluated in quadruple 
precision, while evaluation of the rest continues in double precision. (In some really difficult 
cases, it will be preferable to adopt quadruple precision for the entire integration domain.) 
This strategy has been very successful and many integrals were evaluated very precisely in 
this manner. In fact, in some cases, the achieved numerical precision was such that it led to 
uncovering of errors in some analytic or semi-analytic calculations [[?[]. 

Unfortunately, it is not always easy to detect problems caused by the d-d distortion. This 
was the case with some muon g — 2 diagrams involving sixth-order vacuum-polarization dia- 
grams ||. While Ref. || reported the value — 0.2415(19) (a/vr) 4 the same diagrams evaluated 
using the Pade approximation method gave f|] 

af{Pade) = -0.230362 (5) (^) (15) 

which disagreed with Ref. || by about 6 standard deviations. The reliability of the error 
estimate in || was thus called into question. 

As was discussed above, the most effective way to separate the d-d error evaluated in 
double precision from the statistical uncertainty is to go over to quadruple precision. We 
have therefore repeated the integration of Ref. || entirely in quadruple precision with a 
roughly equal amount of sampling statistics. The new value 



i®{Fig.l) = -0.2285 (18) (^J (16) 

agrees within one standard deviation with the Pade value (|I~5"D but disagrees strongly with 
the old value. 

A closer examination of Ref. || reveals that some of the eight integrals (in particular, the 
diagram e of Fig. 1, evaluated with 60 million sampling points per iteration) show signs of 



oA 4 



suffering from the presence of d-d error. This is now clearly confirmed by the new result (16). 
One way to overcome this problem within double precision calculation is to go over to much 
larger statistics. In order to see whether ten times more sampling points per iteration is 
sufficient for the diagram e, we have evaluated it in both double f] and quadruple^ precisions 
(with 600 million sampling point per iteration and 60 iterations). The two results are in 
good agreement showing that the double precision calculation is reliable now. Encouraged 



2 Evaluated on the Fujitsu VX computer at the Computer Center of Nara Women's University. 
3 Evaluated on the IBM SP2 computer at Cornell Theory Center. 
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by this, we have evaluated the remaining seven integrals in double precision with ten times 
more statistics. The total contribution is found to be 

af{Fig.l) = -0.230 596 (416) (~V (17) 

This is in good agreement with the Pade result fll"5|) . 

We conclude that the problem encountered by the result || was caused solely by insuffi- 
cient statistics. Unfortunately this does not become visible until more extensive calculation 
is done. 

To examine the sensitivity of the Lamb shift calculation to the precision of the arithmetic 
used we evaluated the sixth-order vacuum-polarization contribution to the muonic hydrogen 
in both double and quadruple precision. The results obtained in double and quadruple 
precisions are listed in Table I of Ref. 0. The double precision calculation was carried 
out using 100 million sampling points per iteration while that of quadruple precision was 
obtained for 1 million sampling points per iteration (except for the diagrams a and e which 
employ 2 and 4 times more sampling points, respectively). As is clearly seen from Table I 
of Ref. U , the results in double and quadruple precision are in good agreement and give us 
confidence that this problem does not suffer visibly from the d-d problem. The reason why 
the Lamb shift calculation looks less susceptible to the d-d problem is that, unlike the muon 
g — 2 case, the contribution from large momentum transfer region is strongly suppressed by 
the hydrogenic wave function. 



B. How good is the Pade approximation ? 

We now want to examine whether or not the Pade result a^(Pade) of flip]) agrees within 
its much smaller error bars with the true value of a^(Fig.l). This question is raised because 

the [2/3] or [3/2] Pade approximant P(oS) of (0) does not have complete information on II3 . 
In principle it is possible to answer this question by improving the numerical precision of 
(|T7j) by two orders of magnitude. This is not very practical, however, since it would require 
10 4 -fold increase in the computational effort. Another approach is to go to higher-order 
Pade approximations, taking additional exact properties of njj 1 ' into account P,|10|. 



Let us instead follow an alternative and easier route, namely, examine how well the Pade 
method works by applying it to the cases where exact results are known. 

The first example is the contribution of the two-loop vacuum-polarization to the muon 
g — 2. The exact result is given by ||TTH 



af(a 2 v.p.) = 1.493 671 80 (4) (~) 3 , (18) 

where the uncertainty comes only from the uncertainty in the measured value of m^/m e . 
The corresponding value obtained by numerical integration by VEGAS starting from the 
Kallen-Sabry spectral function is 

a® (K - S) = 1.493 672 7 (40) . (19) 
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The difference between ( |18D and ( |19"D is well within the error bars of (|19"D and will disappear 
as the numerical precision of ( |i~9l) increases. 

Meanwhile, the real and imaginary parts of [3/2] and [2/3] Pade approximation given in 
Appendix A produce 



,(6), 



ay (PadeRe) 



and 



,(6), 



a*-"' (Padelm) 



1.493 677 6 (16) 
1.493 678 1 (16) 
1.493 676 2 (16) 
1.493 675 9 (16) 



1.493 674 2 (38) 
1.493 847 4 (38) 
1.493 737 9 (38) 
1.493 676 5 (38) 



[2/3]A 
[2/3]B 
[3/2]A 
[3/2]B 



[2/3]A 
[2/3]B 
[3/2]A 
[3/2]B 



(20) 



(21) 



respectively, where A and B refer to two different sets of determinations of Pade functions. 



Within the context of this approximation it is not possible to tell which of the results ( 20 ) 
and (|2~1~D are good or bad. Since the [m/n] Pade approximation will get closer and closer 
to the exact function as m and n increase, however, the Pade results based on the real 
and imaginary parts should give identical results in the large m and n limit. This means 
that near-equality of results from the real and imaginary parts may be chosen as a working 
criterion for selecting good Pade approximation. According to this criterion the results 
[2/3]A and [3/2]B are good while [2/3]B and [3/2]A are not. Since this is not a rigorous 
criterion, we may of course choose the worst case as a measure of the uncertainty of the 
Pade approximation. 



The disagreement between better values of (pOj) and (pij ) and the exact value (|i8D is about 
0.000006 (a/n) 3 (~ 0.0004 percent), which falls somewhat outside the estimated errors of 
(0) and A. The difference of about 0.00017 (a/ixf (~ 0.01 percent) between (@) and 



the "worst" value of (|20|) and (|2l|) is much larger, and is well outside its estimated error. 
Based on these observations we may conclude that the near-equality of real and imaginary 
Pade values is a useful working criterion and indicates that the [2/3] and [3/2] Pade methods 
give results, with a high probability, within 0.0004 percent of the exact value. Analogously, 
the Pade approximation value (|i~5"|) of an eighth-order muon g — 2 may not deviate from the 
true value by more than 0.001 percent, even allowing for possible differences between IT2 
and Ilg 1 ', staying within the error bars of ([15]). 

The second example is the Lamb shift caused by the fourth-order vacuum-polarization 
function. The Lamb shift evaluated by VEGAS using the exact vacuum-polarization function 

m is 



AE. 



(p4) 
exact 



0.045 922 738 (57) m r (Zay[ - 



a 



71 



(22) 



while the real and imaginary parts of the [3/2] and [2/3] Pade approximation described in 
Appendix A give 



S 



AP(P4) 
^^PadeRe 



and 



AE 



Padelm 



0.045 923 200 88 
0.045 923 320 16 
0.045 922 988 50 
0.045 922 938 13 



0.045 923 200 97 
0.045 933 403 92 
0.045 926 642 63 
0.045 922 938 22 



m r (Za) 



mJZa) 2 [ — 
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[2/3]A 
[2/3]B 
[3/2]A 
[3/2]B 



[2/3]A 
[2/3]B 
[3/2]A 
[3/2]B 



(23) 



(24) 



respectively^]. Uncertainties of all numerical coefficients in (|23"D and (^4) do not exceed 
0.5 x 10" 9 . 

The [2/3] A and [3/2] B results of fl23|) and <$Z§ are in good agreement with (H), the 
difference being about 0.0007 percent. Meanwhile the [2/3]B value of ([2"4]) deviates from 



by about 0.02 percent. The six standard deviation difference between (22) and "good" 



values of (^f) and (|24j ) is unlikely to be attributable to the fault of numerical integration. 
Instead it must be interpreted as the measure of the approximate nature of the [2/3] and 
[3/2] Pade, which can be remedied only by going to higher-oder Pade approximation. 

Finally one may ask how can one justify the use of the Pade approximation for the 
two-loop vacuum-polarization function as a model for the three-loop case, in view of the 
fact that the three-loop function LT^ has a threshold singularity (~ l/t>) while II 2 does not. 
Our answer is that there is not much analytical difference between P(uj) of two-loop and 
three-loop cases since this singularity is removed explicitly in the construction (|8|) of the 
Pade function P(oj). 



C. Conclusion 

Based on these considerations we conclude: 

(I) The Lamb shift calculation by numerical integration method gives a reliable result 
even if double precision arithmetic is used, provided sufficient sampling of the integrand is 
made in each step of iteration. 

(II) The [2/3] and [3/2] Pade methods for the vacuum-polarization function given in 
Q is a very good approximation over nearly all momentum space region, even though it 
is exact only at a few chosen values of the momentum. It will be possible to narrow the 
gap between Pade and exact results by using more input date and by going to higher-order 
Pade approximations. However, we believe that [2/3] and [3/2] cases are good enough for 
the purpose of confirming the reliability of numerical integration method. 



The results fl23| ) and ( |24| ) are obtained using the numerical integration routine of Maple V. 
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APPENDIX A: PADE APPROXIMATION TO U 2 

The proper two-loop vacuum-polarization function II2 can be obtained from Eq. (57) of 
Ref. m by removing the contribution of (TIi) 2 , where JTi is the one-loop vacuum-polarization 
function. Since n 2 is known explicitly, it can be used to construct a Pade approximant of 
any degree of precision. Our purpose here is to test the approximation to II3 . Thus we 
construct a Pade approximant of IT 2 using analytic information similar to those used in 
constructing described in Sec. [TT]. They are: 
(i) Small momentum behavior of U. 2 (z) given by Q 

„ , , 82 449 2 249916 , 

U 2 (z) = —z + z 2 + z 3 + ■■■ . Al 

v ' 81 675 496125 v ; 



(ii) Leading terms of H 2 (z) for z — > —00: 

3 



n 2 (,) = A-C3~ln(-4,) + i 



ln(-4z) 
4 



+ 0(1). (A2) 



(iii) The threshold behavior of the imaginary part of H 2 (z): 

-Imn 2 (z) = ^ - 25 + 0{5 2 ) , (A3) 
71 ' 4 

where 5 = Jl — l/z. This means that the real part of H 2 (z) has the logarithmic threshold 
singularity of the form 

7T 2 

n 2 (*) = — ln\l-z\ + ■■■ . (A4) 

To take these analytic properties into account we construct a function tl 2 (z) containing 
five parameters ci, c 5 : 

U 2 (z) = IL 2 (z) + (1 - z)G(z){ Cl + c 2 /z) + c 3 + cjz + c 5 F{z), (A5) 

where G(z) is introduced in (|j) and F(z) — In |1 — z\. For z > 1, F(z) has the imaginary 
part 

ImF(z) = -7T . (A6) 

The coefficients Ci, c 2 and C4 are determined by requiring cancellation of all logarithmic 
singularities found in Eqs. (|A2|) and (|A4|) . The remaining C3 and C5 are determined by 



10 



TABLE II. Coefficients of Pade approximants for two- loop vacuum polarization function 112(2). 
We set 6q = 1 fo r the overall normalization. 



Coefficient 


[2/3] A 


[2/3] B 


[3/2]A 


[3/2]B 


a 


0.218 599 301 


0.218 599 301 


0.218 599 301 


0.218 599 301 


a\ 


-0.011 600 062 


-0.097 946 773 


-0.051 634 466 


0.005 017 399 


a 2 


-0.118 517 746 


-0.176 492 097 


-0.190 422 374 


-0.143 212 168 


a 3 








-0.023 342 491 


-0.018 586 700 


bo 


1 


1 


1 


1 


bi 


-0.228 103 636 


-0.623 103 575 


-0.411 244 220 


-0.152 085 728 


b 2 


-0.470 961 931 


-0.667 030 155 


-0.767 838 808 


-0.597 234 583 


b 3 


0.060 545 956 


0.082 509 741 









requiring that LT2 vanishes in the limit 2; — ^ 0_. We find that ^(z) is analytic in the z-plane 
cut along the real axis (1, 00) if we choose 

c a = ^(l-7T 2 ), C2 = i(7-7r 2 ), C4 = l(l + 57r 2 ), 

c 4 = i(-7 + 7r 2 ), c 5 = V (A7) 
The function to be Pade-approximated is chosen as 

' %(z)-%{-oo)\ z=-^—. (A8) 



(1 + ^) 2 V ^ ' zy 7' (1+w 

We construct the Pade approximant P{uj) of the form (|9|) using ( |A1| ), ( |A2| ) and ( |A3[ ) as the 
input data. We obtained two sets of a^s and 6j's for each [2/3] and [3/2] Pade approxima- 
tions, which we denote [2/3] A, [2/3]B, [3/2] A, and [3/2]B. They are listed in Table II. The 
Pade results of Sec. |IIIB| are obtained from this Table. 



As a final check, we compared directly the imaginary parts of the exact Kallen- 
Sabry function and its Pade approximation. The difference of the two functions A = 
Imn 2 (K — S)/tt — Imn 2 (Pade)/7T is examined for the entire momentum range z > 1. The 
fractional deviation A/(Imn 2 (K — S)/n) is less than 0.2 percent in all cases of the Pade 
approximations. The Pade approximation functions are found to oscillate around the exact 



value. Also, by integrating A over the entire range of 5 = Jl — 1/z, < 5 < 1, we found 



/' 

Jo 



0.000 276 • • • [2/3]A 
-0.000 930 • • ■ [2/3]B 



'0 d5 A ~ 1 -0.000 461 • • • [3/2]A (A9) 
0.000 015 ■ ■ • [3/2]B 

Apparently the values of the muon g — 2 and Lamb shift are more accurate than are implied 
by these results. 
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